| Descargar notebook.

5. Compilación de Resultados y Validación

5.1. Cargar las bibliotecas requeridas

El siguiente código carga las bibliotecas requeridas.

import numpy as np                                            # ndarrays for gridded data
import pandas as pd                                           # DataFrames for tabular data

import geostatspy.GSLIB as GSLIB                              # GSLIB utilities, visualization and wrapper
import geostatspy.geostats as geostats                        # GSLIB methods convert to Python
import geostatspy
print('GeostatsPy version: ' + str(geostatspy.__version__))

import os                                                     # set working directory, run executables

from tqdm import tqdm                                         # suppress the status bar
from functools import partialmethod
tqdm.__init__ = partialmethod(tqdm.__init__, disable=True)

ignore_warnings = True                                        # ignore warnings?

import matplotlib.pyplot as plt                               # for plotting
import matplotlib as mpl                                      # custom colorbar
import matplotlib.ticker as mticker                           # custom colorpar ticks
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks
plt.rc('axes', axisbelow=True)                                # plot all grids below the plot elements
if ignore_warnings == True:
    import warnings
    warnings.filterwarnings('ignore')

import os
import matplotlib.image as mpimg
import subprocess
import time
from matplotlib.offsetbox import AnchoredText
import matplotlib.gridspec as gridspec

import plotly.io as pio
import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as py

# trimmed inverted rainbow colormap hsv removing magenta tones
cmap = mpl.colors.ListedColormap(mpl.cm.hsv(np.linspace(0.0, 0.70, 256))).reversed()

GeostatsPy version: 0.0.78

def GSLIB2ndarray_3D(data_file, kcol,nreal, nx, ny, nz):
    if nz > 1 and ny > 1:
        array = np.ndarray(shape = (nreal, nz, ny, nx), dtype=object, order="F")
    elif ny > 1:
        array = np.ndarray(shape=(nreal, ny, nx), dtype=object, order="F")
    else:
        array = np.ndarray(nreal, nx, dtype=object)

    with open(data_file) as f:
        head = [next(f) for _ in range(2)]  # read first two lines
        line2 = head[1].split()
        ncol = int(line2[0])  # get the number of columns

        for icol in range(ncol):  # read over the column names
            head = next(f)
            if icol == kcol:
                col_name = head.split()[0]
        for ineal in range(nreal):
            if nz > 1 and ny > 1:
                for iz in range(nz):
                    for iy in range(ny):
                        for ix in range(nx):
                            head = next(f)
                            array[ineal][nz - 1 - iz][ny - 1 - iy][ix] = head.split()[kcol]
            elif ny > 1:
                for iy in range(ny):
                    for ix in range(0, nx):
                        head = next(f)
                        array[ineal][ny - 1 - iy][ix] = head.split()[kcol]
            else:
                for ix in range(nx):
                    head = next(f)
                    array[ineal][ix] = head.split()[kcol]
    return array, col_name

Nuestro punto de partida sera la base de datos compositada a 6 metros y simplificada pero completa, es decir, sin filtrar aquellas muestras por UG.

## Cargar data simplificada de la sección anterior

import pandas as pd
import numpy as np

# Cargar base de datos en pandas y aproximar a dos decimales
df = pd.read_csv('data.csv', sep=',', encoding='latin1').round(2)
df['UG'] = df['UG'].astype(int)  # Asegurarse de que la columna 'UG' sea de tipo entero

df # mostrar la data

	X	Y	Z	cu_pct	UG
0	472187.34	6925805.69	4211.77	0.02	2
1	472187.74	6925806.53	4205.84	0.01	2
2	472188.09	6925807.38	4199.92	0.01	2
3	472188.47	6925808.28	4194.00	0.01	2
4	472188.88	6925809.22	4188.08	0.01	2
...	...	...	...	...	...
16877	471922.69	6925391.61	4022.52	0.09	4
16878	471924.63	6925394.53	4017.65	0.11	4
16879	471926.59	6925397.45	4012.78	0.18	4
16880	471928.55	6925400.36	4007.92	0.15	4
16881	471930.53	6925403.26	4003.05	0.21	4

16882 rows × 5 columns

# Estadísticas por UG (Number, Max, Min, Mean, Variance)
# Asegúrese de haber ejecutado previamente la celda que carga 'df' (data.csv)
import pandas as pd

# lista de UGs a reportar (ajustar si hay más o diferentes!!!)
ug_list = [2, 3, 4]

stats = {}
for ug in ug_list:
    serie = df.loc[df['UG'] == ug, 'cu_pct']
    stats[f'UG {ug}'] = serie.agg(['count', 'max', 'min', 'mean', 'var'])

# estadisticas globales (todas las muestras)
global_stats = df['cu_pct'].agg(['count', 'max', 'min', 'mean', 'var'])

# montar DataFrame final con el orden solicitado - AJUSTAR SI HAY MÁS UGS!!!
df_stats = pd.DataFrame({
    'UG 2': stats['UG 2'],
    'UG 3': stats['UG 3'],
    'UG 4': stats['UG 4'],
    'Global': global_stats
})

# Renombrar índices a los labels solicitados
df_stats.index = ['Número', 'Máximo', 'Mínimo', 'Media', 'Varianza']

# Mostrar con 2 decimales (ajustar si lo desea)
from IPython.display import display
display(df_stats.round(2))

# Opcional: exportar a csv para revisión fuera del notebook y ajustar decimales
# df_stats.round(6).to_csv('./03_resultados/summary_stats_UGs.csv')

	UG 2	UG 3	UG 4	Global
Número	6203.00	2048.00	8631.00	16882.00
Máximo	2.03	1.64	6.72	6.72
Mínimo	0.00	0.00	0.00	0.00
Media	0.06	0.13	0.25	0.17
Varianza	0.01	0.02	0.05	0.04

5.2. Cargar modelos de UGs y Leyes

El siguiente código carga nuestras estimaciones previas, tanto de UGs como de leyes.

xmin = 471116.24; xmax = 473023.89                  # rango de valores x
ymin = 6924721.57; ymax = 6926164.63                        # rango de valores y
zmin = 2950.02; zmax = 4445.78                              # rango de valores z

xsiz = 15; ysiz = 15; zsiz = 15                             # tamaño de celda
nx = 128; ny = 97; nz = 100                         # número de celdas

tmin = -999; tmax = 999;                                      # límites de recorte de datos

# calculo automático del centro de las capas y cortes para plotear, no modificar

icx = int(nx/2); cx = xmin + icx*xsiz
icy = int(ny/2); cy = ymin + icy*ysiz
icz = int(nz/2); cz = zmin + icz*zsiz

# Cargar modelos de UGs
modelo_ug = GSLIB2ndarray_3D("./03_resultados/modelo_ug.dat ", 0, 1, nx, ny, nz)[0][0].astype(int)

# Cargar modelos de leyes

estCU_UG2 = GSLIB2ndarray_3D("./03_resultados/modelo_leyes_UG2.dat ", 0, 1, nx, ny, nz)[0][0].astype(float)
categoria_UG2 = GSLIB2ndarray_3D("./03_resultados/modelo_categorias_UG2.dat ", 0, 1, nx, ny, nz)[0][0]

estCU_UG3 = GSLIB2ndarray_3D("./03_resultados/modelo_leyes_UG3.dat ", 0, 1, nx, ny, nz)[0][0].astype(float)
categoria_UG3 = GSLIB2ndarray_3D("./03_resultados/modelo_categorias_UG3.dat ", 0, 1, nx, ny, nz)[0][0]

estCU_UG4 = GSLIB2ndarray_3D("./03_resultados/modelo_leyes_UG4.dat ", 0, 1, nx, ny, nz)[0][0].astype(float)
categoria_UG4 = GSLIB2ndarray_3D("./03_resultados/modelo_categorias_UG4.dat ", 0, 1, nx, ny, nz)[0][0]

## Combinar modelos de leyes en un solo modelo final de acuerdo a las UGs
modelo_leyes_final = np.full((nz, ny, nx), np.nan)  # inicializar con NaN
for iz in range(nz):
    for iy in range(ny):
        for ix in range(nx):
            ug_value = modelo_ug[iz, iy, ix]
            if ug_value == 2:
                modelo_leyes_final[iz, iy, ix] = estCU_UG2[iz, iy, ix]
            elif ug_value == 3:
                modelo_leyes_final[iz, iy, ix] = estCU_UG3[iz, iy, ix]
            elif ug_value == 4:
                modelo_leyes_final[iz, iy, ix] = estCU_UG4[iz, iy, ix]

## Combinar modelos de categorías en un solo modelo final de acuerdo a las UGs
modelo_categorias_final = np.full((nz, ny, nx), "Sin_Categoria", dtype=object)  # inicializar con "Sin_Categoria"
for iz in range(nz):
    for iy in range(ny):
        for ix in range(nx):
            ug_value = modelo_ug[iz, iy, ix]
            if ug_value == 2:
                modelo_categorias_final[iz, iy, ix] = categoria_UG2[iz, iy, ix]
            elif ug_value == 3:
                modelo_categorias_final[iz, iy, ix] = categoria_UG3[iz, iy, ix]
            elif ug_value == 4:
                modelo_categorias_final[iz, iy, ix] = categoria_UG4[iz, iy, ix]

continuous = 'cu_pct'

x_coords = np.linspace(xmin + xsiz/2, xmax - xsiz/2, nx)
y_coords = np.linspace(ymin + ysiz/2, ymax - ysiz/2, ny)
z_coords = np.linspace(zmin + zsiz/2, zmax - zsiz/2, nz)

import plotly.io as pio
# pio.renderers.default = "colab"
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode
import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as py

fig = go.Figure(data=[go.Scatter3d(
    x=df['X'],
    y=df['Y'],
    z=df['Z'],
    marker=dict(color=df[continuous],
                colorscale=px.colors.sequential.Rainbow[1:],
                cmin=0.0,
                cmax=df[continuous].quantile(0.95),
                size=2.0, # Esta opcion controla el tamaño de los datos
                colorbar=dict(
                    title='Cu [%]',
                    thickness=20,
                    # Add a border to the colorbar
                    outlinecolor='black',  # Color of the border
                    outlinewidth=2,       # Width of the border
                    bordercolor='white',   # Background border color (if applicable)
                    borderwidth=1         # Background border width
                )
                ),
    mode='markers',
    opacity=1
)])

# otras opciones para controlar el color de los elementos de la figura, borrar hovermode para desactivar etiquetas al pasar el mouse
fig.update_layout(    title_text='Estimación de Cu - UG 4',
scene = dict(
                    xaxis = dict(
                        title='Este',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white",),
                    yaxis = dict(
                        title='Norte',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white"),
                    zaxis = dict(
                        title='Elevación',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white",),
                    ),
                    font_family="Times New Roman",
                    font_color="black", hovermode=False
                  )

# ocultar la leyenda/autonombre "trace 0" junto al colorbar
fig.update_layout(showlegend=False)

# definir cortes en x, y, z

icz = 60  # índice de corte en z (puede ajustarse), va entre 0 y nz-1

# Agregar cortes del modelo de UGs con la misma escala de colores de plotly
x_slice = np.flip(np.flipud(modelo_leyes_final[:, :, icx]), axis=1)
y_slice = np.flip(modelo_leyes_final[:, icy, :], axis=0)
z_slice = np.flip(modelo_leyes_final[nz-icz, :, :], axis=0)

# X-slice (plano constante x)
YY, ZZ = np.meshgrid(y_coords, z_coords)             # shapes (nz, ny)
XX = np.full_like(YY, x_coords[icx])
surf_x = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=x_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=df[continuous].quantile(0.95),
    showscale=False,
    opacity=1,
    name=f'corte X={x_coords[icx]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Y-slice (plano constante y)
XX, ZZ = np.meshgrid(x_coords, z_coords)             # shapes (nz, nx)
YY = np.full_like(XX, y_coords[icy])
surf_y = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=y_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=df[continuous].quantile(0.95),
    showscale=False,
    opacity=1,
    name=f'corte Y={y_coords[icy]:.1f}',
    hoverinfo='skip',
    hovertemplate=None,
)

# Z-slice (plano constante z)
XX, YY = np.meshgrid(x_coords, y_coords)             # shapes (ny, nx)
ZZ = np.full_like(XX, z_coords[icz])
surf_z = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=z_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=df[continuous].quantile(0.95),
    showscale=False,
    opacity=1,
    name=f'corte Z={z_coords[icz]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Añadir superficies al figure (manteniendo la leyenda y colores consistentes)
fig.add_trace(surf_x)
fig.add_trace(surf_y)
fig.add_trace(surf_z)

fig.show()

import plotly.graph_objects as go
import plotly.express as px
from plotly.subplots import make_subplots

ROWS = 3
COLS = 2
# Crear listado de bancos
step_size = 15
max_benches = ROWS*COLS
lims = np.arange(4250,4250-step_size*max_benches,-step_size)

continuous = 'cu_pct'

# Initialize figure with subplots
# Initialize figure with subplots
fig = make_subplots(
    rows=ROWS, cols=COLS, horizontal_spacing=0.2, vertical_spacing=0.1,
    subplot_titles=[f"Cu [%], banco {lims[i]+step_size/2}" for i in range(len(lims))]
)

for i, lim in enumerate(lims):
    data_temp = df.loc[(df['Z'] > lim) & (df['Z'] <= lim + step_size)]
    # Calculate row and column indices for the subplot (1-based indexing)
    row = i // COLS + 1
    col = i % COLS + 1

    # agregar Z-slice del modelo_leyes_final correspondiente al banco (plano constante z)
    z_mid = lim + step_size/2.0
    iz = int(np.abs(z_coords - z_mid).argmin())  # índice de la capa más cercana al centro del banco

    # extraer y orientar la rebanada (debe quedar shape (ny, nx) con y_coords como filas)
    z_slice2d = np.flip(modelo_leyes_final[nz-iz, :, :].astype(float), axis=0)

    # enmascarar valores inválidos (-999) para que sean transparentes en el heatmap
    z_slice2d[z_slice2d <= -998] = np.nan

    heat = go.Heatmap(
        x=x_coords,
        y=y_coords,
        z=z_slice2d,
        colorscale=px.colors.sequential.Rainbow[1:],
        zmin=0.0,
        zmax=df[continuous].quantile(0.95),
        showscale=False,
        opacity=0.7,
        hoverinfo='skip',
        name=f'Z-slice {z_coords[iz]:.1f}'
    )

    # añadir primero el heatmap (fondo) y luego el scatter (muestras) conserva la visibilidad
    fig.add_trace(heat, row=row, col=col)

    fig.add_trace(
        go.Scatter(
            x=data_temp['X'],
            y=data_temp['Y'],
            marker=dict(
                color=data_temp[continuous],
                colorscale=px.colors.sequential.Rainbow[1:],
                cmin=0.0,
                cmax=df[continuous].quantile(0.95),
                size=5.0,
                line=dict(
                    color='black',
                    width=0.5
                ),
            ),
            mode='markers',
            opacity=1,
            showlegend=False
        ),
        row=row, col=col
    )



    # Update xaxis properties
    fig.update_xaxes(title_text="Este", range=[df['X'].quantile(0.015), df['X'].quantile(0.97)], gridcolor='rgba(0,0,0,0.2)', row=row, col=col)

    # Update yaxis properties
    fig.update_yaxes(title_text="Norte", range=[df['Y'].quantile(0.03), df['Y'].quantile(0.99)], gridcolor='rgba(0,0,0,0.2)', row=row, col=col)

fig.update_annotations(font_size=12)
# otras opciones para controlar el color de los elementos de la figura
fig.update_layout(font_family="Times New Roman",font_size=8,
                    font_color="black", plot_bgcolor='white',   autosize = False,
    width = 600,
    height = 900, hovermode=False)


fig.show()

from plotly.subplots import make_subplots

## truncate Z variable in the df axis every x meters

df_swath = df.copy()

dz = 50 # truncar cada 50 metros
df_swath['Z'] = df_swath['Z'].where(df_swath['Z'] % dz == 0, other=df_swath['Z'] - df_swath['Z'] % dz)

## Plot truncated z values against filtered average cut > 0 percentage with a line plot and marginal counter histogram
import plotly.express as px
import plotly.graph_objects as go

# prepare aggregated line data
agg = df_swath[df_swath['cu_pct'] >= 0].groupby('Z').agg(
    mean_cut=('cu_pct', 'mean'),
    count=('cu_pct', 'count')
).reset_index()

# create subplot: top = line, bottom = marginal histogram (shared x)
# enable a secondary y-axis for the bottom subplot so we can plot the bar on the right
fig = make_subplots(rows=2, cols=1, shared_xaxes=True,
                    vertical_spacing=0.05,
                    row_heights=[0.7, 0.3],
                    specs=[[{"type": "xy"}],
                           [{"type": "xy", "secondary_y": True}]])

# line trace (top)
fig.add_trace(
    go.Scatter(
        x=agg['Z'],
        y=agg['mean_cut'],
        mode='lines+markers',
        name='Ley Promedio Cu [%] en muestras',
        line=dict(color='blue'),
        marker=dict(size=6, color='blue')
    ),
    row=1, col=1
)

## add z vs filtered line
# line trace (top)
fig.add_trace(
    go.Scatter(
        x=z_coords,
        y=[np.mean(np.mean(modelo_leyes_final[nz-iz-1, :, :][modelo_leyes_final[nz-iz-1, :, :] > 0])) for iz in range(nz)],
        mode='lines',
        name='Ley Promedio Cu (%) en modelo',
        line=dict(color='red')
    ),
    row=1, col=1
)

# marginal histogram (bottom)
fig.add_trace(
    go.Histogram(
        x=df_swath[df_swath['cu_pct'] >= 0]['Z'],
        nbinsx=30,
        marker=dict(color='blue'),
        opacity=0.3,
        name='Conteo muestras'
    ),
    row=2, col=1
)

# secondary axis: count of valid cells (>0) per model z-layer
valid_counts = [
    int(np.sum((modelo_leyes_final[nz - iz - 1, :, :] > 0) & np.isfinite(modelo_leyes_final[nz - iz - 1, :, :])))
    for iz in range(nz)
]

bar = go.Bar(
    x=z_coords,
    y=valid_counts,
    name='N° celdas válidas (modelo)',
    marker=dict(color='red'),
    opacity=0.3
)
## add bar y-axis on the right, row=2, col=1.

# add the bar to the bottom subplot using the secondary y-axis
fig.add_trace(bar, row=2, col=1, secondary_y=True)
# Update 2nd y-axis properties for the bar chart
fig.update_yaxes(
    title_text='Celdas válidas (modelo)',
    secondary_y=True,
    row=2, col=1
)

# layout and axis labels
fig.update_xaxes(title_text='Elevación (m)', row=2, col=1)
fig.update_yaxes(title_text='Ley Promedio Cu (%)', row=1, col=1)
fig.update_yaxes(title_text='Conteo muestras', row=2, col=1)
fig.update_yaxes(title_text='Celdas válidas (modelo)', row=2, col=1, secondary_y=True)

fig.update_layout(
    title='Ley Promedio de Cu vs Elevación (m)',
    font_family="Times New Roman",
    font_size=12,
    font_color="black",
    plot_bgcolor='white',
    hovermode=False,
    showlegend=True,
    bargap=0.1,
    height=600,
    width=900
)
fig.show()

## esta vez para UG 4

from plotly.subplots import make_subplots

## truncate Z variable in the df axis every x meters

df_swath = df[df['UG'] == 4].copy()

dz = 50 # truncar cada 50 metros
df_swath['Z'] = df_swath['Z'].where(df_swath['Z'] % dz == 0, other=df_swath['Z'] - df_swath['Z'] % dz)

## Plot truncated z values against filtered average cut > 0 percentage with a line plot and marginal counter histogram
import plotly.express as px
import plotly.graph_objects as go

# prepare aggregated line data
agg = df_swath[df_swath['cu_pct'] >= 0].groupby('Z').agg(
    mean_cut=('cu_pct', 'mean'),
    count=('cu_pct', 'count')
).reset_index()

# create subplot: top = line, bottom = marginal histogram (shared x)
# enable a secondary y-axis for the bottom subplot so we can plot the bar on the right
fig = make_subplots(rows=2, cols=1, shared_xaxes=True,
                    vertical_spacing=0.05,
                    row_heights=[0.7, 0.3],
                    specs=[[{"type": "xy"}],
                           [{"type": "xy", "secondary_y": True}]])

# line trace (top)
fig.add_trace(
    go.Scatter(
        x=agg['Z'],
        y=agg['mean_cut'],
        mode='lines+markers',
        name='Ley Promedio Cu [%] en muestras',
        line=dict(color='blue'),
        marker=dict(size=6, color='blue')
    ),
    row=1, col=1
)

modelo_leyes_final_UG4 = np.full((nz, ny, nx), np.nan)  # inicializar con NaN
for iz in range(nz):
    for iy in range(ny):
        for ix in range(nx):
            if modelo_ug[iz, iy, ix] == 4:
                modelo_leyes_final_UG4[iz, iy, ix] = estCU_UG4[iz, iy, ix]


## add z vs filtered line
# line trace (top)
fig.add_trace(
    go.Scatter(
        x=z_coords,
        y=[np.mean(np.mean(modelo_leyes_final_UG4[nz-iz-1, :, :][modelo_leyes_final_UG4[nz-iz-1, :, :] > 0])) for iz in range(nz)],
        mode='lines',
        name='Ley Promedio Cu (%) en modelo',
        line=dict(color='red')
    ),
    row=1, col=1
)

# marginal histogram (bottom)
fig.add_trace(
    go.Histogram(
        x=df_swath[df_swath['cu_pct'] >= 0]['Z'],
        nbinsx=30,
        marker=dict(color='blue'),
        opacity=0.3,
        name='Conteo muestras'
    ),
    row=2, col=1
)

# secondary axis: count of valid cells (>0) per model z-layer
valid_counts = [
    int(np.sum((modelo_leyes_final_UG4[nz - iz - 1, :, :] > 0) & np.isfinite(modelo_leyes_final_UG4[nz - iz - 1, :, :])))
    for iz in range(nz)
]

bar = go.Bar(
    x=z_coords,
    y=valid_counts,
    name='N° celdas válidas (modelo)',
    marker=dict(color='red'),
    opacity=0.3
)
## add bar y-axis on the right, row=2, col=1.

# add the bar to the bottom subplot using the secondary y-axis
fig.add_trace(bar, row=2, col=1, secondary_y=True)
# Update 2nd y-axis properties for the bar chart
fig.update_yaxes(
    title_text='Celdas válidas (modelo)',
    secondary_y=True,
    row=2, col=1
)

# layout and axis labels
fig.update_xaxes(title_text='Elevación (m)', row=2, col=1)
fig.update_yaxes(title_text='Ley Promedio Cu (%)', row=1, col=1)
fig.update_yaxes(title_text='Conteo muestras', row=2, col=1)
fig.update_yaxes(title_text='Celdas válidas (modelo)', row=2, col=1, secondary_y=True)

fig.update_layout(
    title='Ley Promedio de Cu vs Elevación (m) - UG 4',
    font_family="Times New Roman",
    font_size=12,
    font_color="black",
    plot_bgcolor='white',
    hovermode=False,
    showlegend=True,
    bargap=0.1,
    height=600,
    width=900
)
fig.show()

from IPython.display import display

# entregar estadísticas de modelo_leyes_final según UG (Number, Max, Min, Mean, Variance)

ugs = [2, 3, 4]

results = {}
for ug in ugs:
    mask = (modelo_ug == ug) & (modelo_leyes_final > 0)
    vals = modelo_leyes_final[mask]
    cnt = vals.size
    if cnt > 0:
        vals = vals.astype(float)
        results[f'UG {ug}'] = {
            'Número': int(cnt),
            'Máximo': float(np.nanmax(vals)),
            'Mínimo': float(np.nanmin(vals)),
            'Media': float(np.nanmean(vals)),
            'Varianza': float(np.nanvar(vals))
        }
    else:
        results[f'UG {ug}'] = {'Número': 0, 'Máximo': np.nan, 'Mínimo': np.nan, 'Media': np.nan, 'Varianza': np.nan}

# Estadísticas globales (todas las celdas válidas)
mask_all = (modelo_leyes_final > 0) & np.isfinite(modelo_leyes_final)
vals_all = modelo_leyes_final[mask_all]
cnt_all = vals_all.size
if cnt_all > 0:
    results['Global'] = {
        'Número': int(cnt_all),
        'Máximo': float(np.nanmax(vals_all)),
        'Mínimo': float(np.nanmin(vals_all)),
        'Media': float(np.nanmean(vals_all)),
        'Varianza': float(np.nanvar(vals_all))
    }
else:
    results['Global'] = {'Número': 0, 'Máximo': np.nan, 'Mínimo': np.nan, 'Media': np.nan, 'Varianza': np.nan}

# Montar DataFrame y mostrar
df_stats_model = pd.DataFrame(results).T[['Número', 'Máximo', 'Mínimo', 'Media', 'Varianza']]
display(df_stats_model.round(2))

	Número	Máximo	Mínimo	Media	Varianza
UG 2	73198.0	0.70	0.00	0.06	0.00
UG 3	13488.0	0.74	0.01	0.13	0.01
UG 4	124567.0	2.48	0.01	0.23	0.01
Global	211253.0	2.48	0.00	0.16	0.02

# Ensure renderer is appropriate for your environment (notebook / vscode)
pio.renderers.default = 'notebook'

# --- Input assumptions (these must exist in the environment) ---
# df: pandas DataFrame with columns 'X','Y','Z'
# categoria_UG4: numpy array with shape (nz, ny, nx) containing strings: 'Inferido','Indicado','Medido'
# x_coords, y_coords, z_coords: 1D coordinate arrays for grid cell centers (lengths nx, ny, nz)
# icx, icy, icz: integer indices for the central slices to display
# nz, ny, nx: grid dimensions

# Map categories to numeric codes
cat_to_num = {'Inferido': 0, 'Indicado': 1, 'Medido': 2}
cat_num = np.vectorize(lambda x: cat_to_num.get(x, np.nan))(modelo_categorias_final)
cat_num = cat_num.astype(float)  # allow NaNs

# Create slices (same convention used in the notebook)
x_slice = np.flip(np.flipud(cat_num[:, :, icx]), axis=1)   # X-slice (nz, ny)
y_slice = np.flip(cat_num[:, icy, :], axis=0)             # Y-slice (nz, nx)
z_slice = np.flip(cat_num[nz-icz, :, :], axis=0)          # Z-slice (ny, nx)

# Discrete colors: Inferido (red), Indicado (yellow), Medido (green)
# We use normalized stops [0.0, 0.5, 1.0] so numeric codes 0,1,2 map exactly to the three colors
colorscale = [
    [0.0, '#d73027'],   # Inferido -> 0
    [0.5, "#fff700"],   # Indicado -> 1 (normalized to 0.5)
    [1.0, '#1a9850']    # Medido -> 2 (normalized to 1.0)
]

# Build surface for X-slice
YY, ZZ = np.meshgrid(y_coords, z_coords)             # shapes (nz, ny)
XX = np.full_like(YY, x_coords[icx])
surf_x = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=x_slice.astype(float),
    colorscale=colorscale,
    cmin=0, cmax=2,
    showscale=False,
    opacity=1,
    name=f'corte X={x_coords[icx]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Build surface for Y-slice
XX, ZZ = np.meshgrid(x_coords, z_coords)             # shapes (nz, nx)
YY = np.full_like(XX, y_coords[icy])
surf_y = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=y_slice.astype(float),
    colorscale=colorscale,
    cmin=0, cmax=2,
    showscale=False,
    opacity=1,
    name=f'corte Y={y_coords[icy]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Build surface for Z-slice (we enable a colorbar here with categorical ticks)
XX, YY = np.meshgrid(x_coords, y_coords)             # shapes (ny, nx)
ZZ = np.full_like(XX, z_coords[icz])
surf_z = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=z_slice.astype(float),
    colorscale=colorscale,
    cmin=0, cmax=2,
    showscale=True,
    colorbar=dict(
        title='Categoría',
        tickvals=[0, 1, 2],
        ticktext=['Inferido', 'Indicado', 'Medido'],
        thickness=20
    ),
    opacity=1,
    name=f'corte Z={z_coords[icz]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Add sample points as neutral black markers for reference
scatter_samples = go.Scatter3d(
    x=df['X'], y=df['Y'], z=df['Z'],
    mode='markers',
    marker=dict(color='black', size=2, opacity=0.8),
    name='Muestras'
)

# Compose figure
fig = go.Figure(data=[scatter_samples, surf_x, surf_y, surf_z])
fig.update_layout(
    title='Categoría (Inferido=rojo, Indicado=amarillo, Medido=verde)',
    scene=dict(
        xaxis=dict(title='Este'),
        yaxis=dict(title='Norte'),
        zaxis=dict(title='Elevación')
    ),
    width=900,
    height=720,
    font=dict(family='Times New Roman', size=12)
)

fig.show()

# Entregar para cada UG y cada categoría las estadísticas (Número, Max, Min, Mean, Variance) de modelo_leyes_final
from IPython.display import display

ugs = [2, 3, 4]
categories = ['Inferido', 'Indicado', 'Medido']
results = {}
for ug in ugs:
    for cat in categories:
        mask = (modelo_ug == ug) & (modelo_categorias_final == cat) & (modelo_leyes_final > 0)
        vals = modelo_leyes_final[mask]
        cnt = vals.size
        key = f'UG {ug} - {cat}'
        if cnt > 0:
            vals = vals.astype(float)
            results[key] = {
                'Número': int(cnt),
                'Máximo': float(np.nanmax(vals)),
                'Mínimo': float(np.nanmin(vals)),
                'Media': float(np.nanmean(vals)),
                'Varianza': float(np.nanvar(vals))
            }
        else:
            results[key] = {'Número': 0, 'Máximo': np.nan, 'Mínimo': np.nan, 'Media': np.nan, 'Varianza': np.nan}
# Montar DataFrame y mostrar
df_stats_cat = pd.DataFrame(results).T[['Número', 'Máximo', 'Mínimo', 'Media', 'Varianza']]
display(df_stats_cat.round(2))

	Número	Máximo	Mínimo	Media	Varianza
UG 2 - Inferido	1430.0	0.43	0.00	0.04	0.00
UG 2 - Indicado	8147.0	0.53	0.00	0.05	0.00
UG 2 - Medido	63621.0	0.70	0.00	0.06	0.00
UG 3 - Inferido	295.0	0.43	0.03	0.08	0.00
UG 3 - Indicado	1024.0	0.44	0.01	0.10	0.00
UG 3 - Medido	12169.0	0.74	0.01	0.13	0.01
UG 4 - Inferido	2195.0	0.55	0.03	0.20	0.01
UG 4 - Indicado	19573.0	0.68	0.03	0.22	0.01
UG 4 - Medido	102799.0	2.48	0.01	0.24	0.01

# Entregar para cada UG y cada categorías global = med+ind+inf y demostrado = med + ind las estadísticas (Número, Max, Min, Mean, Variance) de modelo_leyes_final
from IPython.display import display

ugs = [2, 3, 4]
categories = [['Inferido', 'Indicado', 'Medido'], ['Indicado', 'Medido']]
category_labels = ['Global', 'Demostrado']
results = {}

for ug in ugs:
    for cat, label in zip(categories, category_labels):
        mask = (modelo_ug == ug) & np.isin(modelo_categorias_final, cat) & (modelo_leyes_final > 0) & np.isfinite(modelo_leyes_final)
        vals = modelo_leyes_final[mask]
        cnt = vals.size
        key = f'UG {ug} - {label}'
        if cnt > 0:
            vals = vals.astype(float)
            results[key] = {
                'Número': int(cnt),
                'Máximo': float(np.nanmax(vals)),
                'Mínimo': float(np.nanmin(vals)),
                'Media': float(np.nanmean(vals)),
                'Varianza': float(np.nanvar(vals))
            }
        else:
            results[key] = {'Número': 0, 'Máximo': np.nan, 'Mínimo': np.nan, 'Media': np.nan, 'Varianza': np.nan}
# Montar DataFrame y mostrar
df_stats_cat = pd.DataFrame(results).T[['Número', 'Máximo', 'Mínimo', 'Media', 'Varianza']]
display(df_stats_cat.round(2))

	Número	Máximo	Mínimo	Media	Varianza
UG 2 - Global	73198.0	0.70	0.00	0.06	0.00
UG 2 - Demostrado	71768.0	0.70	0.00	0.06	0.00
UG 3 - Global	13488.0	0.74	0.01	0.13	0.01
UG 3 - Demostrado	13193.0	0.74	0.01	0.13	0.01
UG 4 - Global	124567.0	2.48	0.01	0.23	0.01
UG 4 - Demostrado	122372.0	2.48	0.01	0.23	0.01

## Curva Tonelaje (izq.) y ley promedio (der.) in secondary axis para el modelo final de leyes overlayed in the same plot
import plotly.graph_objects as go
# Prepare data
cutoff_values = np.arange(0, df['cu_pct'].quantile(0.96), 0.01)
mtonnage = []
average_grade = []
for cutoff in cutoff_values:
    mask = (modelo_leyes_final >= cutoff) & np.isfinite(modelo_leyes_final)
    ## puede crear otra mascara mask para filtrar por categoría o UG si lo desea
    mask = np.logical_and(mask, np.logical_or(modelo_categorias_final == 'Medido', modelo_categorias_final == 'Indicado', modelo_categorias_final == 'Inferido'))  # ejemplo: solo categoría Medido
    ## mask = np.logical_and(mask, modelo_ug == 4)  # ejemplo: solo UG 4. COMBINAR AMBOS FILTROS SI ES NECESARIO

    total_tonnage = np.sum(mask) * (xsiz * ysiz * zsiz * 2.5)  # assuming 2.5 t/m³ (2.5 g/cm³)
    total_mtonnage = total_tonnage / 1e6  # convert to million tonnes
    mtonnage.append(round(total_mtonnage, 2))
    if total_tonnage > 0:
        avg_grade = np.nanmean(modelo_leyes_final[mask])
    else:
        avg_grade = 0
    average_grade.append(round(avg_grade, 2))

# Create figure with secondary y-axis
fig = go.Figure()
# Tonnage trace
fig.add_trace(
    go.Scatter(
        x=cutoff_values,
        y=mtonnage,
        name='Tonelaje',
        line=dict(color='blue'),
        yaxis='y1'
    )
)
# Average grade trace
fig.add_trace(
    go.Scatter(
        x=cutoff_values,
        y=average_grade,
        name='Ley Promedio (%)',
        line=dict(color='red'),
        yaxis='y2'
    )
)

# Update layout for dual y-axes
fig.update_layout(
    title='Curva Tonelaje y Ley Promedio vs Corte de Ley (Modelo Final de Leyes)',
    xaxis=dict(title='Corte de Ley (%)',gridcolor='rgba(0,0,0,0.2)'),
    yaxis=dict(
        title='Tonelaje (MTon)',
        tickfont=dict(color='blue',
                      ),
        side='left',
        gridcolor='rgba(0,0,0,0.2)'
    ),
    yaxis2=dict(
        title='Ley Promedio (%)',
        tickfont=dict(color='red'),
        tickformat=',.f',
        overlaying='y',
        side='right'
    ),
    legend=dict(x=0.1, y=0.9),
    font_family="Times New Roman",
    font_size=12,
    font_color="black",
    plot_bgcolor='white',
    width=800,
    height=500
)
fig.show()